#!/usr/bin/perl
# 2013/05/24 Make graph of distance vs association frequnecy

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 6 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [map file] -o [output file] -c [combination number file]\n";
}

my %opt;
getopts("i:o:c:", \%opt);
my $FILE_map = $opt{i};
my $FILE_out = $opt{o};
my $FILE_combi = $opt{c};

my $Resolution = 5000;

#---------------------------------------
# read combination information
#---------------------------------------
my %COMBINATION;
my $fh_combi = IO::File->new($FILE_combi) or die "cannot open $FILE_combi: $!";
while($_ = $fh_combi->getline()){
	s/\r?\n//;
	my ($distance, $number) = split /\t/;
	$COMBINATION{$distance} = $number;
}
$fh_combi->close();


#---------------------------------------
# collect information
#---------------------------------------
my %data;
my $fh_map = IO::File->new($FILE_map) or die "cannot open $FILE_map: $!";
while($_ = $fh_map->getline()){
	s/\r?\n//;
	my ($id, $chr1, $loc1, $direction1, $mapQ1, $hinID1, $hinLoc1, $chr2, $loc2, $direction2, $mapQ2, $hinID2, $hinLoc2, $type, $num, $locations) = split /\t/;

	if($chr1 ne $chr2){
		next;
	}

	if($type ne 'U:U'){
		next;
	}

	my $distance= abs($loc1 - $loc2);
	my $category = int($distance / $Resolution) * $Resolution;
	$data{$category}++;
}
$fh_map->close();


#---------------------------------------
# output distance vs contact probability
#---------------------------------------
my $fh_out = IO::File->new($FILE_out, 'w') or die "cannot write $FILE_out: $!";
$fh_out->print("distance\tobserved\ttotal_combination\tcontactProbability\n");
foreach my $distance(sort{$a <=> $b} keys %data){
	my $observed = $data{$distance};
	my $combination = exists $COMBINATION{$distance} ? $COMBINATION{$distance} : 0;
	my $probability = $combination == 0 ? 0 : $observed / $combination;
	$fh_out->print("$distance\t$observed\t$combination\t$probability\n");
}
$fh_out->close();

